/*
 * DxFile.cpp
 *
 *  Created on: Nov 9, 2011
 *      Author: yi
 */

#include "mbuild/DxFile.h"

DxFile::DxFile() {
	g=NULL;
}

DxFile::~DxFile() {
	delete g;
}
DxFile::DxFile(GridSpace *g)
{
	load(g);
}

void DxFile::load(GridSpace *g)
{
	this->g=new GridSpace;
	g->copyto(this->g);
}

DxFile::DxFile(char *fname)
{
	g=new GridSpace;
	DxFile::read_dx(fname,g);

}

void DxFile::read(char *fname)
{
	g= new GridSpace;
	DxFile::read_dx(fname,g);
}



void DxFile::write(char *fname)
{
	if (g == NULL) {
		throw EmptyDxFile();
	}
	DxFile::write_dx(fname,g);
}



void DxFile::read_dx(char *fname, GridSpace *readto)
{
	FILE* in;
	char str[MAX_LINE_SIZE],s[MAX_LINE_SIZE];
	double tmp,x0,y0,z0,dx,dy,dz;
	unsigned long ncel_x,ncel_y,ncel_z,ncel,tmpi;
	in = fopen(fname,"r");
	if (in == NULL) {
		throw FileNotOpen();
	}

	/* First read the header */
	do {
		fgets(str,MAX_LINE_SIZE,in);
	} while (strchr(str,'#')==str);
	sscanf(str, "%6s %1s %5s %13s %6s %lu %lu %lu \n", s,s,s,s,s,&ncel_x,&ncel_y,&ncel_z);
	readto->set_dimension(ncel_x,ncel_y,ncel_z);
	fgets(str,MAX_LINE_SIZE,in);
	sscanf(str, "%6s %lf %lf %lf\n",s, &x0, &y0, &z0);
	readto->set_origin(x0,y0,z0);
	fscanf(in, "%5s %lf %lf %lf \n",s, &dx, &tmp, &tmp);
	fscanf(in, "%5s %lf %lf %lf \n",s, &tmp, &dy, &tmp);
	fscanf(in, "%5s %lf %lf %lf \n",s, &tmp, &tmp, &dz);
	readto->set_grid_size(dx,dy,dz);
	fgets(s,MAX_LINE_SIZE,in);
	fscanf(in, "%6s %lu %5s %5s %4s %6s %4s %lu %5s %lu %4s %7s \n",s,&tmpi,s,s,s,s,s,&tmpi,s,&ncel,s,s);
	readto->initialize();

	double* d=readto->values();
	double x,y,z;
	/* Read in the rest of data */
	int i,last,nline,offset;
	last = ncel%3; // the number of datapoint in last line
	nline =(ncel-last)/3;    /* total lines  */
	offset =0;

	for (i=1; i <= nline; ++i)
	{
	fscanf(in,"%lf %lf %lf \n", &x, &y, &z);
	d[offset]=x;
	d[offset+1]=y;
	d[offset+2]=z;
	offset+=3;
	}

	if (last == 1)
	fscanf(in,"%lf \n", &d[offset]);    /* reading in the last line */
	else if (last == 2)
	fscanf(in,"%lf %lf \n", &d[offset], &d[offset+1]);
	fclose(in);
}





void DxFile::ave(char *outfilename)
{
	//Write the average values on Z axis to a output file.
	FILE* out;
	double x0,y0,z0,dx,dy,dz,x,y,z;
	out = fopen(outfilename,"w");
	if (out == NULL) {
		throw FileNotOpen();
	}
	if (g == NULL ) {
		throw EmptyDxFile();
	}
	unsigned long i,ncel,nx,ny,nz,iz;
	ncel=g->get_dimension(nx,ny,nz);
	double* ave= new double[nz];
	double* std= new double[nz];
	double* data=g->values();
	g->get_grid_size(dx,dy,dz);
	g->get_origin(x0,y0,z0);
	for(i=0;i<nz;i++) {ave[i]=0.0; std[i]=0.0;}
	for(i=0;i< ncel;i++) {
		iz=i%nz;
		ave[iz]+=data[i]/(nx*ny);
		std[iz]+=data[i]*data[i]/(nx*ny);
	}
	fprintf(out,"## Z ave  std\n");
	for(i=0;i<nz;i++) {
		z=i*dz+z0;
		std[i]-=ave[i]*ave[i];
		fprintf(out,"%7.4f %7.4f %7.4f\n",z,ave[i],sqrt(std[i]));
	}
	delete [] ave;
	delete [] std;
	fclose(out);
}

void DxFile::rave(char *outfilename)
{
	FILE* out;
	double x0,y0,z0,dx,dy,dz,x,y,z,r;
	out = fopen(outfilename,"w");
	if (out == NULL) {
		throw FileNotOpen();
	}
	unsigned long i,j,ncel,nx,ny,nz,iz,ir;
	ncel=g->get_dimension(nx,ny,nz);
	double* data=g->values();
	g->get_grid_size(dx,dy,dz);
	g->get_origin(x0,y0,z0);

	if(dx!= dy) throw XYNotEqual();
	double dr=sqrt(dx*dy);
	int nr=int(dx*nx/(dr*2));
	double** ave= new double*[nr];
	double** std= new double*[nr];
	unsigned long** nave= new unsigned long*[nr];
	for(i=0;i <nr;i++) {
		ave[i] = new double[nz];
		std[i] = new double[nz];
		nave[i] = new unsigned long[nz];
	}
	for( i=0;i<nr;i++)
	{
		for( j=0;j<nz;j++)
		{
			ave[i][j]=0;
			nave[i][j]=0;
			std[i][j]=0;
	    }
	}
	for( i=0;i<ncel;i++)
	{
		x=(i/(nz*ny))*dx+x0;
		y=((i/nz)%ny)*dy+y0;
		z=(i%nz)*dz+z0;
		r=sqrt(x*x+y*y);
		ir=floor(r/dr);
		iz=i%nz;
		if (ir<nr) {
			ave[ir][iz]+=data[i];
	        std[ir][iz]+= data[i]*data[i];
	        nave[ir][iz]++;
		}
	}
	fprintf(out,"## R Z ave std\n");
	  for( i=0;i<nr;i++)
	   {
	      for(j=0;j<nz;j++)
	      {
	    	  if(nave[i][j]>0) {
	    		  ave[i][j]/=nave[i][j];
	    		  std[i][j]/=nave[i][j];
	    		  std[i][j]-=ave[i][j]*ave[i][j];
	    		  std[i][j]=sqrt(std[i][j]);
	    		  r=i*dr;
	    		  z=z0+j*dz;
	    		  fprintf(out,"%7.4f %7.4f %7.4f %7.4f\n",r,z,ave[i][j],std[i][j]);
	          }
	        }
	   }

	  fclose(out);
	  for(i=0;i <nr;i++) {
			delete [] ave[i];
			delete [] std[i];
			delete [] nave[i];
		}
	  delete [] ave;
	  delete [] std;
	  delete [] nave;
}


void DxFile::write_dx(char *fname, GridSpace *writefrom)
{
	/* Write gridspace to file*/

	 FILE* out;
	 double x0,y0,z0,dx,dy,dz;
	 unsigned long ncel_x,ncel_y,ncel_z,ncel;
	 out = fopen(fname,"w");


	 if (out == NULL) {
		 throw FileNotOpen();
	}
	 ncel=writefrom->get_dimension(ncel_x,ncel_y,ncel_z);

	 //If GridSpace is empty, do nothing.
	 if (ncel <=0) {
	 		throw EmptyDxFile();
	 }
	 /* write the header part of DX*/
	 fprintf(out, "object 1 class gridpositions counts %lu %lu %lu \n",ncel_x,ncel_y,ncel_z);
	 writefrom->get_origin(x0,y0,z0);
	 fprintf(out, "origin %12.6E %12.6E %12.6E \n", x0, y0, z0);
	 writefrom->get_grid_size(dx,dy,dz);
	 fprintf(out, "delta %12.6E %12.6E %12.6E \n", dx,0.000000E+00,0.000000E+00);
	 fprintf(out, "delta %12.6E %12.6E %12.6E \n", 0.000000E+00,dy,0.000000E+00);
	 fprintf(out, "delta %12.6E %12.6E %12.6E \n", 0.000000E+00,0.000000E+00,dz);
	 fprintf(out, "object 2 class gridconnections counts %lu %lu %lu\n", ncel_x, ncel_y, ncel_z);
	 fprintf(out, "object 3 class array type double rank 0 items %lu data follows\n", ncel_x*ncel_y*ncel_z);

	/* Write the data points */
	double* map=(writefrom->values());
	unsigned long i,last,nline,offset;
	last = ncel%3; //number of data point in last line
	nline =(ncel-last)/3;    /* total lines */
	offset =0;
	for (i=1; i <= nline  ; ++i)
	{
	  fprintf(out, "%12.6E %12.6E %12.6E \n", cleandata(map[offset]), cleandata(map[offset+1]), cleandata(map[offset+2]));
	offset+=3;
	}

	if (last == 1)
	  fprintf(out,"%12.6E \n", cleandata(map[offset]));    /* saving in the last line */
	else if (last == 2)
	  fprintf(out,"%12.6E %12.6E \n", cleandata(map[offset]),cleandata( map[offset+1]));
	fclose(out);
}



inline double DxFile::cleandata(double number)
{
	/*  clean the data to prevent scrambled number */
	  if (number<1.000000e-99 && number>-1.000000e-99 ) return 0.000000e+00;
	  else return number;
}



